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Abstract 



Large scale molecular dynamics (MD) simulations on two-phase immisci- 
ble flows show that associated with the moving contact line, there is a very 
large 1/x partial-slip region where x denotes the distance from the contact 
line. This power-law partial-slip region is verified in large-scale adaptive con- 
tinuum simulations based on a local, continuum hydrodynamic formulation, 
which has proved successful in reproducing MD results at the nanoscale. Both 
MD and continuum simulations indicate the existence of a universal slip pro- 
file in the Stokes- flow regime, well described by v shp (x)/V w = 1/(1 + x/al s ), 
where v shp is the slip velocity, V w the speed of moving wall, l s the slip length, 
and a is a numerical constant. Implications for the contact-line dissipation 
are discussed. 
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Boundary condition that specifies the flow of fluid over a solid surface is a cornerstone 
of hydrodynamics. The no-slip condition, i.e., zero relative velocity between the fluid and 
solid at the interface, has been the paradigm in most of the hydrodynamics literature [1]. 
In molecular dynamics (MD) simulations, however, a small amount of relative slip between 
the fluid and the solid surface is generally detected at high flow rate [2]. Such slip can 
be accounted for by the Navier boundary condition (NBC), whereby the slip velocity is 
proportional to the tangential viscous stress [2,3]. The proportionality constant between 
the slip velocity v shp and the shear rate is denoted the slip length l s , usually ranging from 
one to a few nanometers (from MD simulations). As the amount of slip is extremely small 
in subsonic flow rates, the NBC is practically indistinguishable from the no-slip condition 
in most situations. In contrast, for immiscible flows the MD simulations have shown near- 
complete slip in the vicinity of the moving contact line (MCL) , defined as the intersection of 
fluid-fluid interface with the solid wall [4-6]. An intriguing question ensues: In a mesoscopic 
or macroscopic system, what is the slip profile which consistently interpolates between the 
near complete slip at the MCL and the no-slip boundary condition that must hold at regions 
far away [7,8]? Recent evidences have shown the slip profile obtained from nanoscale MD 
simulations to be accountable by the generalized Navier boundary condition (GNBC), in 
which the slip velocity is proportional to the total tangential stress - - the sum of the 
viscous stress and the uncompensated Young stress; the latter arises from the deviation of 
the fluid-fluid interface from its static configuration [6]. Here we show through large-scale 
MD simulations and continuum hydrodynamic calculations, that there exists a power-law 
partial-slip region extending to hundreds of nanometers or even more, contrary to the usual 
expectation of a nanometer-scale slip region in the vicinity of the MCL. The existence of 
this large partial-slip region modifies the conventional picture that significant NBC slipping 
occurs only under high flow/shear rate. Instead, the power-law slip region is associated with 
the universal slip profile of the MCL, even at low flow rates. 

MD simulations have been carried out for increasingly larger systems of immiscible 
Couette flow (Fig. 1). Two immiscible fluids were confined between two parallel walls 
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in the xy plane, with the fluid-solid boundaries defined by z = 0, H . Periodic bound- 
ary conditions were imposed along the x and y directions. Interaction between fluid 
molecules separated by a distance r was modeled by a modified Lennard- Jones (LJ) po- 
tential Uff = 4e (a/r) 12 — 5ff (cr/r) 6 , where <5// = 1 for like molecules and 5ff = —1 
for molecules of different species. The average number density for the fluids was set at 
p = 0.81/cr 3 . The temperature was controlled at 2.8e/ ks- Each wall was constructed by 
two [001] planes of an fee lattice, with each wall molecule attached to a lattice site by a 
harmonic spring. The mass of the wall molecule was set equal to that of the fluid molecule 
m. The number density of the wall was set at p w = 1.86/a 3 . The wall-fluid interaction 
was modeled by another LJ potential U w f with the energy and range parameters given by 
e w f = 1.16e and a w f = 1.04a, and S w f for specifying the wetting property of the fluid, taken 
to be 1 in our simulations. The Couette flow was generated by moving the top and bottom 
walls at a constant speed V w in the ±x directions, respectively. In most of our simulations, 



the shearing speed was V w = 0.05^/e/m, the sample dimension along y was 6.8a, the wall 
separation along z varied from H = 6.8a to 68a, and the sample dimension along x was set 
to be long enough so that the uniform single-phase shear flow was recovered far away from 
the MCL. Steady-state interfacial and velocity profiles were obtained from time average over 



5 x 10 5 r where r is the atomic time scale yma 2 /e. 

The tangential slip velocity profiles next to the wall, i.e., the slip profiles, are shown in 
the inset to Fig. 2. While there is clearly a small core region, on the order of a few l s , 
where the slip profiles display sharp decay, a much more gentle variation of the slip profiles 
becomes apparent as the system size H increases. In order to quantify the nature of the 
gentle variation, we plot in Fig. 2 the same data in the log-log scale. The dashed line has 
a slope —1, indicating the 1/x behavior to be indeed realized in MD simulations. For our 
finite-sized systems, there is always a plateau in the slip velocity in each of the single-phase 
flows, also observable in MD simulations with a value given by v'^ hp = 2V w l s /(H + 2l s ), 
which acts as an outer cutoff on the 1/x profile. This v'q Hp expression is simply derived 
from the Navier-Stokes equation for uniform shear flow and the NBC. From our largest MD 



simulation, the v shp oc 1/x behavior is seen to extend to ~ 50a (or ~ 25Z S ) . Hence as 
H — > oo and t>Q hp approaches (no-slip), the power-law region can be very large indeed. 
A large 1/x partial slip region is significant, because the outer cutoff length scale directly 
determines the integrated effects, such as the total steady-state dissipation. While in the 
past the 1/x stress variation away from the MCL has been known [9], to our knowledge the 
fact that the partial slip also exhibits the same spatial dependence has not been previously 
seen [10], even though the validity of the Navier boundary condition at high shear stress has 
been verified [2,3]. 

Since MD simulations reach size and accuracy (e.g., for V w / ' \Jefm <C 1 ) limits quickly, 
a continuum hydrodynamic formulation is necessary for realistic simulations. In particular, 
continuum simulations are necessary for low flow rate immiscible flows, where MD simula- 
tions are known to be very resource intensive. Combining the GNBC with the Cahn-Hilliard 
(CH) hydrodynamic formulation of two-phase flow [11,12], we have obtained a continuum 
hydrodynamic model [6], suitable for the calculation of much larger immiscible flows (than 
MD simulations) that are accurate to the molecular scale. In the same notations of Ref. [6], 
the continuum model is formulated as follows. The two coupled equations of motion are the 
Navier-Stokes (NS) equation (with the addition of the capillary force density) and the CH 
convection-diffusion equation for the composition field 0(r) = (p 2 — Pi)/(p2 + Pi) (where p\ 
and p 2 are the local number densities for the two fluid species): 

= -Vp + V • <j v + pV0 + mpg ext , (la) 

^ + v-V0 = MV 2 /i, (lb) 

together with the incompressibility condition V • v = 0. Here mp is the average fluid mass 
density, p is the pressure, cr v is the Newtonian viscous stress tensor, /iV0 is the capillary 
force density with p = 8F/8(j) being the chemical potential defined from the CH free energy 
functional F [13], mpg ex t is the external body force density (for Poiseuille flows), and M is 
the phenomeno logical mobility coefficient. The boundary conditions at the solid surface are 
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v n = 0, d n fi — (n denotes the outward surface normal), the continuum form of the GNBC: 

Pvfp = -r)d n v x + L{4>)d x 4>, (lc) 
and the relaxational equation for surface 0: 

^ + vV0 = -rL(0), (id) 

Here L(<f>) = Kd n 4> + d^ w f(4>) / d<fi with 7«,/(0) being the wall- fluid interfacial free energy 
density, L(<j))d x <j) is the uncompensated Young stress, and T is a (positive) phenomenological 
parameter. 

The continuum results shown in Fig. 2 were calculated on a uniform mesh, using the 
same set of material parameters and M, T values corresponding to the same local properties 
in all the five MD simulations. The overall agreement is excellent. Such agreement is possible 
because the GNBC does not impose an artificial cutoff on the slip region. Below we extend 
the MCL simulations, through continuum hydrodynamics, to lower flow rates and much 
larger systems. 

The continuum simulation for macroscopic immiscible flows is a challenging task. Meth- 
ods based on a fixed uniform mesh would break down because it cannot afford to simulate 
macroscopic systems with molecular resolution near the MCL. We have employed for this 
problem the adaptive method based on iterative grid redistribution [14]. The computa- 
tional mesh is redistributed according to the behavior of the continuum solution so that fine 
molecular resolution is achieved in the interfacial region and near the MCL, while elsewhere 
a much coarser mesh is used to save computational cost. The mesh distribution is controlled 
by a monitor function (see [14]), and the redistribution procedure is done repeatedly as the 
solution evolves to its steady state. A semi-implicit time stepping scheme is also used to 
speed up the approach to steady state. For Couette flow, since the interface and contact 
line are confined in a region that is narrow in the x direction but extended in the entire z 
direction, we used a variable, adaptive grid in the x direction and kept a uniform grid in the 
z direction. This relatively simpler grid structure (compared to a general two dimensional 
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variable grid) greatly simplifies the space discretization as well as the matrix structure in 
the implicit time discretization. 

Figure 3 shows the continuum results for three large systems (H on the order of hundreds 
of a) with small V w (~ 0.001-^e/m, well beyond the accuracy of our MD simulations). In 
all three cases, the capillary force was verified to be important only in the interfacial region. 
In fact it decays to zero (exponentially) within a few a. However, the pressure gradients 
and viscous forces shows a much slower variation. They are balanced outside the interfacial 
region, indicating the flow to be governed by the Stokes equation. This is expected, because 
the Reynolds number mpV w H/r] pa 0.6 for p pa 0.8/a 3 , V w = 0.005y / e/m, H pa 300<r, and 
rj ~ 2.0y / em/(T 2 . In Fig. 3 the slip profiles, plotted on the log-log scale, clearly show the 
1/x behavior extending from \x\ pa 6l s to \x\ pa 270l s . The inset to Fig. 3 shows the scaled 
tangential velocity profiles at the solid surface, from which the existence of universal slip 
profile is evident. Physically, when H 3> l s , the regime of Stokes flow is governed by only 
one velocity scale V w and one length scale l s . Thus universality becomes evident in terms 
of v x /V w plotted as a function of x/l s . (The inset of Fig. 2 has already shown part of 
the universality, that the core slip profile is independent of the system size H.) We give a 
heuristic account of the universal slip profile as follows. Away from the MCL, the viscous 
shear stress is given by —ar)v x (x)/\x\, where a is a constant ~ 1, rj the viscosity, and v x (x) the 
local tangential velocity. The NBC implies v s J; ip (x) = —l s av x (x)/\x\. Since v l x hp (x) = v x + V w , 
combining the two equations yields v s J ip (x)/V w = 1/(1 + \x\/al s ) [15]. This relation, with 
a pa 2.14 for best fit, agrees with the continuum slip profiles extremely well, as seen in the 
inset to Fig. 3. 

In the Couette geometry, external work is supplied to maintain the constant speed V w 
of the moving wall. The rate of work is given by the integral of the local tangential force 
lj 1 rj\v shp \ times the wall speed V w [6], i.e., 

J (C^l^l) V w dx = J)VIA 

per unit length along y, where A is a numerical constant. In the limit of V w — > and 
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l s /H — * 0, A is independent of V w and l s but depends on the outer cutoff of the l/x profile. 
As the external work done in the steady state must be fully dissipated, the total dissipation 
rate is equal to the rate of external work. 

Slip profiles obtained from both MD and continuum simulations show that the l/x 
partial-slip region starts from x c « 6l s (v shp xs 0.2QV W ). The outer cutoff for the partial- 
slip region, denoted by R, is determined by the overall size of the system. For the total 
dissipation, the contributions of the core region and the partial-slip region may be quantified 
by the dimensionless integrations 
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for the l/x region. Assuming l s ~ 1 nm (a few cr's) [2], for x c = 10 nm and R — 1 /im we 
have V ~ 9.2. While for x c = 10 nm and R = 1 mm we have V ~ 24, i.e., the power-law 
region can contribute significantly more to the total dissipation than the core region. 
The dissipation component that occurs at the fluid-solid interface can be evaluated as 



J (C^l^l 2 ) dx = r^lT 



per unit length along y, where 
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In the core region, the interfacial component of the dissipation is obtained by letting R = x c , 
orI« 1.6, i.e., about 70% of the total interfacial dissipation (X = 2.14 for R/l s — > oo). 

Recently, there has been considerable interest in the transport at micro- and nanoscales. 
The lower limit for R can reach submicrometer or even shorter length scales. On the other 
hand, there have been increasing evidences for large slip length realized in various fluid-solid 
interfaces [16-18]. Slip length l s as large as 1 /im has been reported [17,18]. The results 
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in this paper indicate that fluid-solid interfacial dissipation is an important contribution to 
the total dissipation if large slip length occurs in a small system. While asymptotic analysis 
has shown that at large distances from the MCL, the flow field is not sensitive to the slip 
boundary condition [19], yet the (macroscopic) asymptotic region may not be attained given 
the small system size and/or the large slip length. In this regard, a continuum hydrodynamic 
formulation of the contact-line motion is necessary for realistic simulations of fluid dynamics 
at micro- and nanoscales, as done in the present case. 

Partial support from HKUST's EHIA funding and Hong Kong RGC grants HKUST 
6176/99P, 6 143/0 IP, and 604803 is hereby acknowledged. 
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FIG. 1. Segment of the MD simulation sample for an immiscible Couette flow with a 90° static 
contact angle. The small empty and solid circles indicate the instantaneous molecular positions 
of the two fluids projected onto the xz plane. The horizontally aligned black/gray circles denote 
the wall molecules. Between the two fluids, the large hollow circles represent the time-averaged 
interface profile, defined by pi = P2 (0 = 0). The solid line is the interface profile calculated from 
the continuum hydrodynamic model. 
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FIG. 2. Log- log plot of the slip profiles showing the 1/x behavior. Here v x /V w + 1 is the 
scaled slip velocity at the lower fluid-solid interface z = 0, and x/a measures the distance from 
the MCL in units of a. The wall is moving at — V w , hence v x /V w = means complete slip and 
v x /V w = — 1 means no slip. The v x profiles were obtained for five symmetric cases of identical 
wall-fluid interactions for the two fluids (both 5 w f = 1 and thus a 90° static contact angle). The 
five cases shown here used different values for H but the same value for V w (= 0.05i/e/m) and also 
the same parameters for densities and interactions. The symbols represent the MD results and the 
solid lines represent the continuum hydrodynamics results, obtained for H = 6.8a (black circles 
and line), H = 13.6cr (red squares and line), H = 21.2a (green diamonds and line), H = 54. 4a 
(blue up-triangles and line), H = 68a (orange down-triangles and line). The dashed line has the 
slope of —1, indicating that the 1/x behavior is approached for increasingly larger H. For H = 68a, 
the 1/x behavior extends from \x\ ~ 12a ~ 6l s to 50<r ~ 25l s , where l s was measured to be 2a. 
Inset: The scaled tangential velocity v x /V w at z = is plotted as a function of x/a. 
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FIG. 3. Log-log plot of the slip profiles showing the 1/x behavior. Here v x /V w + 1 is the scaled 
slip velocity at the lower fluid-solid interface z = 0, and x/l s measures the distance from the MCL 
in unit of l s . The wall is moving at — V w , hence v x /V w = means complete slip and v x /V w = — 1 
means no slip. The black solid line denotes the case of H = 326<r, V w = 0.005 y^e/m and l s = 1.24a, 
the red dashed line denotes the case of H = 326a, V w = 0.0025-\/e/m. and l s = 1.24cr, and the 
blue dotted line denotes the case of H = 326cr, V w = 0.0025ve/m and l s = 0.62a. The green 
dot-dashed line has the slope of —1, indicating a power-law region much wider than that in Fig. 
2. Inset: Universal slip profile. The scaled tangential velocities v x /V w at z = for all three cases 
are plotted as a function of the scaled coordinate x/l s . It is seen that the slip profiles show a 
partial-slip region as large as hundreds of l s . The relation v x /V w = 1/(1 + |x|/2.14Z s ) — 1 is also 
plotted by the orange dot-dashed line, showing an extremely good fit. 
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